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Abstract 

We consider an elliptic PDE in two variables. As one parameter approaches zero, 
this PDE collapses to a parabolic one, that is forward parabolic in a part of the do- 
main and backward parabolic in the remainder. Such problems arise naturally in 
various stochastic models, such as fluid models for data-handling systems and Markov- 
modulated queues. We employ singular perturbation methods to study the problem 
for small values of the parameter. 

1 Introduction 

We consider the following boundary value problem for F = F(x,£) 

DF xx + (a-OF x + F^ + ^F) ( = 0; x > 0, -oo < f < oo (1.1) 
DF t (p,Z) + (a-Z)F(0,£) = 0, -oo < f < oo (1.2) 

F(x,£)dxd£ = l. (1.3) 




-oo 



Here a and D are positive constants, and F is a bivariate probability density. This problem 
arose as an asymptotic limit of a Markov-modulated queueing model, where the input process 
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is generated by N sources that turn "on" and "off" at exponential waiting times. When on, a 
source generates a Poisson arrival stream to the queue. The joint distribution of the number 
of on sources and the queue length satisfies a complicated system of difference equations, 
that as N — > oo may be approximated by the problem p.lj) - (jl.3jl . The variable x is related 
to the queue length and £ corresponds to a scaled, centered measure of the number of on 
sources. An empty queue has x — 0, and £ = means that the number of active sources 
equals its mean value. 

The diffusion coefficient D in (jl.lj) measures variability effects in the service time distri- 
bution, while a measures the difference between average output and input rates to the queue. 
The condition a > guarantees stability and the existence of a steady state distribution. 
Note that D — > means that the queue becomes a deterministic, or "fluid", process. A 
detailed derivation of (|l.l|) - (|1.3|l can be found in [T]. There we also showed that the solution 
can be reduced to either a Fredholm integral equation of the second kind (for F(0, £)), or to 
finding a single eigenvector of an infinite matrix, whose elements are expressed in terms of 
Laguerre functions. 

The complexity of the exact solution to (|l.ljl - ([1.3j) suggests that an asymptotic analysis 
may be fruitful. In lj we considered the limit D — > oo and discussed numerical methods 
that work well for D large and D = 0(1). The purpose of this note is to analyze the opposite 
limit, namely D — > + . This is a very singular limit, as we show in what follows. Also, the 
numerical methods involve truncating an infinite matrix, say to an M x M matrix. For 
D — > + the eigenvector we seek decays on the scale (see [1, section 4.2]) 0(D~ 2 ). Thus 
for say D = .1, M must be of the order of about 500 before an accurate result is obtained. 
The basic matrix is not sparse; it does become diagonally dominant for D — > oo, but not for 
D -> 0+ 

Let us briefly discuss the problem with D = 0. Denoting by J-(x, £) the solution to (jl.lj) 
with D = 0, we see that the elliptic PDE degenerates into a parabolic one, that is forward 
parabolic in the range £ > a and backward parabolic for £ < a. The boundary condition (jl.2j) 
becomes .F(0, £) = and can only be applied in the range £ > a, where the PDE is forward 
parabolic. The study of backward/forward parabolic problems dates as far back as 1914 (see 
Gerrey j2J) and particular problems of the type here are analyzed in [3,4,5]. In j3] we give an 
explicit expression for !F(x, £), subject to the condition .F(oo,£) = (27r)~ 1 / 2 e~'> I 2 (replacing 
(jl.3jl ). The function T can also be interpreted as an asymptotic limit of a probability 
distribution for a discrete stochastic model, that was formulated and analyzed in jE]. Here 
we are interested in how the solution of the elliptic PDE approaches that of the parabolic 
one. For the latter the values of ^"(0, £) for £ < a are unknown and must be computed as 
a part of the solution. Probabilistically, this corresponds to boundary mass along this part 
of the boundary. For any D > 0, (jl.3j) shows that F is a proper density function. Now the 
"no-flux" boundary condition (jl.2j) applies for all £, and no boundary mass develops. We 
also note that by integrating (jl.lj) over x G (0, oo) and using (jl.2j) . the marginal distribution 
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Figure 1: A sketch of the deterministic trajectories for a 
J °° F(x,£)dx satisfies an elementary ODE that is readily solved to give 



F(x,0<fe = -7=e -c/2 . (1.4) 
V 2n 
o 

Here we also used the normalization ()1.3|) . Since F is a density in x while is a distribution 
in x, we expect to compare F to JF^. for x > 0. 

To get a qualitative picture of the solution, we plot in Figure 1 the deterministic approx- 
imation to (jl.lj) . 

This corresponds to neglecting diffusion in both the x- and ^-variables, and is the phase- 
plane flow x = £ — a, £ = — £. We plot a few trajectories for x > 0, and note that once 
they hit a; = (necessarily with £ < a), they stay at a; = and flow to the equilibrium at 
£ = 0. The figure also indicates the boundary and corner layers that arise in the analysis for 
D -> 0. 



2 Singular Perturbation Analysis 

We shall analyze (|1.1J1 - (|1.3J) for D — > + in the three ranges (i) x > 0, — oo < £ < oo (or 
x ps 0, £ > a); (ii) 1 » 0, ( < a; and (iii) i « 0, ^ « a. 
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2.1 Outer Solution 



We assume x > and expand F(x,£) as F(x,£) = F (x,£) + DFi(x,£) + 0(D 2 ) to find that 

te-a)F Q>x = F Qi K + (£F )t. (2.1) 

We note that when £ > a, ()2.1|) is a diffusion equation with x taking the place of time, and 
when £ < a it is a diffusion equation run in reverse time. As is well-known [2,3,5], we need to 
impose boundary conditions for (|2.1|) on the "incoming" half of the boundary, where x = 
and £ > a. We can derive the appropriate "half boundary condition" for (|2.1j) by considering 
a boundary layer about x — 0. 

We thus consider the boundary layer where x = 0(D). Setting x = Du and F(x, £) = 
G(it,0 ~ J D I/0 G (u,0, we find that 

G ,uu + (a- OG ,u = 0; u > 0, -oo < £ < oo (2.2) 
G 0jU (0, + (o - £)Gq(0, £) = 0, -oo < £ < oo (2.3) 

and hence Go is proportional to exp [(£ — a)u] . For £ > a this grows exponentially as u — > oo 
and cannot match to the outer solution Fq(x, £) as x — > + , for any choice of i/q- We thus 
conclude that no boundary layer develops along £ > a, where (|2.1|) is forward parabolic. 
Thus Fq must satisfy F (0, £) = for £ > a, which is the limit of (jl.2|) as D — > + . 

The explicit solution to (|2.1jl . subject to the above "half boundary condition," is given 
in There we solved (j2.1)l using a Laplace transform in the x- variable. Due to the for- 
ward/backward nature of the PDE, the spectrum has both positive and negative eigenvalues. 
Insisting that the solution decay as x — > oo, the positive eigenvalues must be absent. This, 
together with the half boundary condition for £ > a is sufficient to determine the solution. 

Noting that Fq should correspond to J- x , we give two different representations of this 
solution. The first is the contour integral 

F (x,0 = K±-.\ ^J^D a e + e^ + 29)A(ey x dB, (2.4) 

9 9(9 + a) 



A(9) = -^exp 
+ a 



~9C{\) +\l9(9 + a) 



m=l 



S n = - [a + Va 2 + An ] , n = 0,1,2,.... 



Here D p (-) is the parabolic cylinder function of order p, 7 is Euler's constant and £(•) is the 
Riemann zeta function. The Bromwich contour Br in (J2.4)) is any vertical contour in the 
half-plane Re(6 l ) > —a, where the integrand is analytic. Note that the only singularities of 
the integrand are simple poles at 9 — —S n , n > 0. If we close the contour to the left, then 
we can convert the integral to the sum of the residues at the poles. This leads to a second 
representation of the solution in terms of the problem's eigenfunctions: 
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K 22 °nt 



-S n X 



-e/2 



n=0 



2n 



a2 _ n/ 2 e _ 5 2 

a - S n 
n > 1, 



exp 



'r 1 

n ^2^ + 2™ 7 



v/2 

oo 

n 

J m=l 



(2.5) 



exp 



m 2m 



Co 



ae exp 



«C(;) 



n 

m=l 



5n 



exp 



is the n th Hermite 



The product in the expression for c n omits the term m = n and H n 
polynomial. 

The constant K is a normalization constant that will ultimately be determined from 
(jl.HJ) . We cannot determine it at this stage since we shall show that there is also an 0(1) 
probability mass in a boundary layer near x = (with £ < a) , that contributes to (jl.3j) . 



2.2 Boundary Layer 

While no boundary layer develops for £ > a, there is one for £ — a < and x = 0(D). We 
again scale x = Du, and set 

F(x, = GK = iGo(u, + Gx(«, + 0(D). (2.6) 
Then G will satisfy the problem (|2.2j) and (|2.3|) . whose general solution is 

G (u,0=9o(0e^- a)u , (2.7) 

and this decays exponentially as u — > oo for £ < a. We will determine <?o(£) shortly, and also 
show that the expansion (J2.6|) must include the term of order 0(D~ l ) in order to match to 
the outer expansion. 

Using ()2.6|) in (jl.l|) we find that the correction term G\ satisfies 

Gi,uu + ( a ~ O^i," = — ^o,k — (^Oo)g (2.8) 
= -e (? " a)u K + (£ + 2u)g' + (u 2 + < + 1)<? ] 

and (|E2D yields 

Gi,„(0,0 + (o-0Gi(0, = 0. (2.9) 
We write the solution to ()2.8|) as 

G 1 (u,0 = e^ a)u G(u,0 + h(0- (2.10) 

Then the BC Q yields _ 

G u (0,O + (a-OMO = 0. (2.11) 
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From (|2SJ and (EHTH) we obtain 

G uu + (f - a)G u = - K + (f O )' + w(2^ + £# ) + u 2 g ] 

so that 



m 3 „,^u 2 



G(u, = + + 7(0« + 5(e) (2-12) 

o 2 



where 

2 e 



£ - a £ - a 



<7o(0, (2-13) 



7(0 = -^W+(^)l+^^[Vo + eflb]-^^sfld 



and <5(0 is an arbitrary function. 

The asymptotic matching of the inner and outer solutions requires that as u — > oo (|2.6|) 
agrees with F as x — > + , and since £ < a this implies that 

MO = F (0,0 = tf-^ / ^4^W(£ + 20M(0)d0. (2.14) 

27TZ y V27T 
Br 

Since G u (0,O = 7(0, (EH yields 

*<>-«->»«>- J « 215 > 

where the last equality follows from ()2.13|) . We thus set go(0 — ( a ~ 0#(0 to find that 

9"(0 + (^)'(O + (a - 0^o(0,0 = 0. (2.16) 
Also, the leading term in the boundary layer becomes 

G(u,0~^(a-0^~ a)u 9(0, (2.17) 

By integrating (|2.1jl from x = to x = oo and comparing the result to (j2.15j) . we conclude 
that P(0 = J °° F (x, i)dx + g(£) satisfies P"(0 + (C^)'(O = for £ < a, so we have 

^(0 + / F (x, Odx = ^=e~ e /\ £ < a (2.18) 

v 27T 



for some constant L. But the left side of (I2.18j) is the leading term in the expansion of the 
left side of (II .4j) for D —>■ + , hence L — 1. Note also that the normalization condition in 
(I1.3J1 implies to leading order that 

a oo oo 

Jg(Z)dt+J J F (x,Z)dxdS=l. (2.19) 

— oo — oo 



It remains only to determine K. For £ > a there is no boundary layer contribution for 
x = 0(D) to CCH) and thus 



/ F (x,t)dx = tf-L / ^=ZW(£ + m^r-M = 4r e ~ e/2 ( 2 - 20 ) 
J 2m J Bl i y/2'K v y/2n 

o 

where we restrict —a < Re(#) < on Br' (this allows us to perform the ^-integration using 
(12. 4j) by simply replacing e by 1/0) . For £ > a we can close the integrating contour in the 
right half-plane, picking up the residue from the pole at 9 = 0. We conclude that K — 1, as 
A(Q) = 1 and £> (0 = e" f2/4 . 

To summarize, we have obtained the outer solution as ()2.4j) with K — 1 and the boundary 
layer as (J2.17)) , with (J2.18|) and L = 1 . The correction term in the boundary layer is given by 
G\ in (J2.6|) . with (|2.1U|) - (|2.14|) . The function 5(£) can ultimately be found by higher order 
matchings; this would require we compute the correction term F\ in the outer solution. The 
analysis shows that for D small and x > (or x ~ with £ > a) the shape of the density 
for the elliptic problem is, to leading order, the same as that for the parabolic one. However 
for x small and £ < a a boundary layer develops for the elliptic problem, and this takes the 
place of the boundary mass present in the parabolic problem. 

2.3 Corner Layer 

Form (j2.17|) we see that the leading term vanishes as £ f a, indicating another non-uniformity 
in the asymptotics. We thus study how the boundary layer disappears as £ increases through 
the critical value a, by introducing the scaling 

£ - a = D 1/A z, x = D 3/4 X (2.21) 

with 

F(x, = H(X, z) ~ D U H (X, z) (2.22) 
where v is a constant to be determined. From (jl.lj) and (jl.2|) we find that Hq satisfies 

H ,xx - zH 0tX + H 0tZZ = 0; X > 0, -oo < z < oo (2.23) 
#q,x(0,2) -£#0(0,3) = 0, -oo<z<oo. (2.24) 

We note that u = D~ l l 4 X and (£ — a)u = zX. Note also that going from the u to X scales 
represents a "thickening" of the boundary layer as £ | a. 

We seek solutions of ()2.23|) in the separable form Hq = e rX TL(z;T), where r G C is a 
separation constant. We find that 7i satisfies the Airy equation 

H zz = (ZT - T 2 )H, 

so that H will be proportional to Ai(r 1 / 3 (z — r)J (we expect Hq to decay as z — > 00, 
corresponding to the disappearance of the boundary layer). We thus argue that a general 
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solution to (|2.23|) will have the form 

H (X, z) = -^J e TX Ai{r^z - r^)f{r)dr (2.25) 

c 

for some function / and contour C in the r-plane. If the contour C can be chosen in the 
range Re(r) < then we integrate (|2.23|) and use (|2.24|) and (|2.25j) to obtain 

so that the integral must be a linear function of z, and / satisfies a Fredholm integral equation 
of the first kind. 

While we have not been able to determine f(r) explicitly, we can show that (|2.25|) can 
asymptotically match with the outer solution in (|2.4j) . For X and z large and positive we 
expect to have to scale r in ()2.25|) to be small. If /(r) ~ fiT Ul as r — > the right side of 
(|2.25|) may be approximated by 

— f fiT Ul e TX Ai(T 1,z z)dT. (2.27) 
2m J 

c 

We next examine ()2.4j) for x — + and £ — > a. For s^Owe must consider 9 — > oo. In ^E] 
we evaluated the infinite product in ()2.4j) as 9 — > oo, with the result 

A(9) ~ ^e C3 e- a2 / 4 v^e-( e2+ ^ lose e e2 / 2 . (2.28) 

Here C3 = 6*3(0) is a constant. With the scaling b — > 00, F — > 00 and F — 2^ = 0(6~ 1 / 6 ), 
the parabolic cylinder function -Db(F) may be approximated by an Airy function: 

D b (Y) ~ 2 fe / 2 r b 1/e Ai(b l ^(Y -2Vb)). (2.29) 

Applying QTZfy with b = a9 + 9 2 and F = f + 26* we have 

(F - 2v / 6) 1/6 = [f + 29 - 2V9 2 + a9][a9 + 9 2 } 1/6 ~ (f - a)^/ 3 

and then 

+ 20) ~ ^ 1 /3 e (V2)(^+a e )log(^+ae) e -(l/ 2 )(^) Ai (^ _ ^1/3) q 3q) 

This holds in the limit -> oo with £ - a = O(0- 1 / 3 ). Combining (j2~2l^|) and (j2~37H) in (jUD 
and recalling that = 1 leads to 

FoM „ ±/ Ae , e -,,«!) i9 ,, 31) 



Br 



ae C3 e- a2 / 4 



2ni / rV6 



Br 
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By comparing D u x (|2.27|) with (|2.3H1 we see that matching is possible with 



-~, * = -~ h = ae c *e-« 2 /\ 
8 o 



(2.32) 



Thus the density in this corner region will be 0(D 5 / 8 ). This means that the total mass in 
the corner is (noting that dxd£ = DdzdX) 0(D 3 / 8 ). 

Finally we investigate the matching between the corner and boundary layers. Noting 
that represents the mass in the boundary layer for a given £ < a, we expect that g — > 
as £ | a - Let us assume that 



0(0 ~ 02(o-O 



i/ 2 



6 -»• a. 



Then since a — £ = — D 1 l^z the boundary layer becomes 0\jD 1 (a — 0" 2+1 ] = 0[D^ U2 3 ^ 4 ] 
on the z scale. This can match to the corner layer if (z/ 2 — 3)/4 = —5/8, i.e., z/ 2 = 1/2. Now 
suppose we let z — ► —00 and X — > in ()2.24j) . scaling r = z + cj(— z) -1 / 3 . This yields 



D-^ S H (X, 



D -5/8 e zX { 



-1/3. 



2m 



^/(-,)V3 /(z+a;( _ z)1/3) (2 ^ 3) 



x ^{[-^-cu(-^)- 1 / 3 ] 1/3 ^(-^) 1 / 3 }^. 
If f(z) ~ / 3 (— 2)^ as z — > —00 the left side of ()2.33|) should approach 



D 8 < : V ./';( : 
This can match to ()2.16|) if 



IAj-l/3 



2tt? 



Az(u;)oL; 



(2.34) 



^3 
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02 = /3 



2vri 



Ai(o;)da; 



This shows that the matchings between the corner layer and other two expansions are possible 
if /( T ) = 0{t~ x / & ) as t -> and /(r) = 0(|r| n / 6 ) as r -> -00. 



3 Conclusion 

Our analysis showed that the density F(x,£) is 0(1) for x > 0, and is large (0(-D -1 )) in 
the boundary layer where x = 0(D) and £ < a. It was also shown to be large (0(D~ 5 ^ 8 )) 
in the corner layer where (x, ~ (0, a), with the precise scaling given by (j2.21j) . The mass 
in the boundary layer and outer region are comparable, while that in the corner layer is 
asymptotically small. We compare this to the solution of the forward/backward problem in 
jlj, where there is non-zero boundary mass along x = and £ < a. 
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We are presently investigating the problem (|l.l|) - ()1.3j) in the limit a — > oo, with D 
fixed. This corresponds, after some rescaling, to having small diffusion in both the x and 
£ directions. Some preliminary results show that a geometrical optics type expansion is 
fruitful, and that the problem yields interesting asymptotic structures, such as a caustic 
boundary and interior cusped caustics. 
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